Paquetes

# instalamos los paquetes necesarios
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("feedeR")) install.packages("feedeR")
if(!require("fs")) install.packages("fs")
if(!require("lubridate")) install.packages("lubridate")
if(!require("fs")) install.packages("fs")
if(!require("tmap")) install.packages("tmap")
if(!require("classInt")) install.packages("classInt")

# cargamos los paquetes
library(feedeR)
library(sf) 
library(fs)
library(tidyverse)
library(lubridate)
library(classInt)
library(tmap)

Enlaces de descarga

La primera url nos dará acceso a un listado de provincias, sedes territoriales (no siempre coinciden con la provincia), con nuevos enlaces RSS los cuales incluyen los enlaces finales de descarga para cada municipio. En este caso, descargaremos el edificado de Valencia. Los datos del Catastro se actualizan cada seis meses.

url <- "http://www.catastro.minhap.es/INSPIRE/buildings/ES.SDGC.bu.atom.xml"

# importamos los RSS con enlaces de provincias
prov_enlaces <- feed.extract(url)
str(prov_enlaces) #estructura es lista

## List of 4
##  $ title  : chr "Download service of Buildings. Territorial Office"
##  $ link   : chr "http://www.catastro.minhap.es/INSPIRE/buildings/ES.SDGC.BU.atom.xml"
##  $ updated: POSIXct[1:1], format: "2019-10-26"
##  $ items  :'data.frame': 52 obs. of  4 variables:
##   ..$ title: chr [1:52] "Territorial office 02 Albacete" "Territorial office 03 Alicante" "Territorial office 04 Almería" "Territorial office 05 Avila" ...
##   ..$ date : POSIXct[1:52], format: "2019-10-26" ...
##   ..$ link : chr [1:52] "http://www.catastro.minhap.es/INSPIRE/buildings/02/ES.SDGC.bu.atom_02.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/03/ES.SDGC.bu.atom_03.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/04/ES.SDGC.bu.atom_04.xml" "http://www.catastro.minhap.es/INSPIRE/buildings/05/ES.SDGC.bu.atom_05.xml" ...
##   ..$ hash : chr [1:52] "d21ebb7975e59937" "bdba5e149f09e9d8" "03bcbcc7c5be2e17" "8a154202dd778143" ...

# extraemos la tabla con los enlaces
prov_enlaces_tab <- as_tibble(prov_enlaces$items)
prov_enlaces_tab

## # A tibble: 52 x 4
##    title          date                link                         hash    
##    <chr>          <dttm>              <chr>                        <chr>   
##  1 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ d21ebb7~
##  2 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ bdba5e1~
##  3 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ 03bcbcc~
##  4 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ 8a15420~
##  5 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ 7d3fd37~
##  6 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ 9c08741~
##  7 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ ff722b1~
##  8 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ b431aa6~
##  9 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ f79c656~
## 10 Territorial o~ 2019-10-26 00:00:00 http://www.catastro.minhap.~ d702a6a~
## # ... with 42 more rows

Accedemos y descargamos los datos de Valencia. Para encontrar el enlace final de descarga usamos la función filter() del paquete dplyr buscando el nombre de la sede territorial y posteriormente el nombre del municipio en mayúsculas con la función str_detect() de stringr. La función pull() nos permite extraer una columna de un data.frame.

# filtramos la provincia y obtenemos la url RSS
val_atom <- filter(prov_enlaces_tab, str_detect(title, "Valencia")) %>% pull(link)

# importamos la RSS
val_enlaces <- feed.extract(val_atom)

# obtenemos la tabla con los enlaces de descarga
val_enlaces_tab <- val_enlaces$items

# filtramos la tabla con el nombre de la ciudad
val_link <- filter(val_enlaces_tab, str_detect(title, "VALENCIA")) %>% pull(link)
val_link

## [1] "http://www.catastro.minhap.es/INSPIRE/Buildings/46/46900-VALENCIA/A.ES.SDGC.BU.46900.zip"

Descarga de datos

La descarga se realiza con la función download.file() que únicamente tiene dos argumentos principales, el enlace de descarga y la ruta con el nombre del archivo. En este caso hacemos uso de la función tempfile(), que nos es útil para crear archivos temporales, es decir, archivos que únicamente existen en la memoría RAM por un tiempo determinado. El archivo que descargamos tiene extensión *.zip, por lo que debemos descomprimirlo con otra función (unzip()), que requiere el nombre del archivo y el nombre de la carpeta donde lo queremos descomprimir. Por último, la función URLencode() codifica una dirección URL que contiene caracteres especiales.

# creamos un archivo temporal 
temp <- tempfile()

# descargamos los datos
download.file(URLencode(val_link), temp)

# descomprimimos a una carpeta llamda buildings
unzip(temp, exdir = "buildings")

Importar los datos

Para importar los datos utilizamos la función dir_ls() del paquete fs, que nos permite obtener los archivos y carpetas de una ruta concreta al mismo tiempo que filtramos por un patrón de texto (regexp: expresión regular). Aplicamos la función st_read() del paquete sf al archivo espacial de formato Geography Markup Language (GML).

# obtenemos la ruta con el archivo
file_val <- dir_ls("buildings", regexp = "building.gml")

# importamos los datos
buildings_val <- st_read(file_val)

## Reading layer `Building' from data source `C:\Users\xeo19\Documents\GitHub\blogR_update\content\post\es\2019-11-01-visualizar-crecimiento-urbano\buildings\A.ES.SDGC.BU.46900.building.gml' using driver `GML'
## Simple feature collection with 36296 features and 24 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 720608 ymin: 4351287 xmax: 734982.5 ymax: 4382906
## epsg (SRID):    25830
## proj4string:    +proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Preparación de los datos

Únicamente convertimos la columna de la edad del edificio (beginning) en clase Date. La columna de la fecha contiene algunas fechas en formato --01-01 lo que no corresponde a ninguna fecha reconocible. Por eso, reemplazamos el primer - por 0000.

# 
buildings_val <- mutate(buildings_val, 
               beginning = str_replace(beginning, "^-", "0000") %>% 
                            ymd_hms() %>% as_date()
               )

## Warning: 4 failed to parse.

Gráfico de distribución

Antes de crear el mapa de la edad del edificado, lo que reflejará el crecimiento urbano, haremos un gráfico de distribución de la fecha de construcción de los edificios. Podremos identificar claramente períodos de expansión urbana. Usaremos el paquete ggplot2 con la geometría de geom_density() para este objetivo. La función font_add_google() del paquete sysfonts nos permite descargar e incluir familias de fuentes desde Google.

#descarga de fuente
sysfonts::font_add_google("Montserrat", "Montserrat")

#usar showtext para fuentes
showtext::showtext_auto() 

#limitamos al periodo posterior a 1750
filter(buildings_val, beginning >= "1750-01-01") %>%
 ggplot(aes(beginning)) + 
    geom_density(fill = "#2166ac", alpha = 0.7) +
  scale_x_date(date_breaks = "20 year", 
               date_labels = "%Y") +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 12, family = "Montserrat"),
        title = element_text(family = "Montserrat"),
        axis.text = element_text(family = "Montserrat")) +
  labs(y = "",x = "", title = "Evolución del desarrollo urbano")

Buffer de 2,5 km de Valencia

Para poder visualizar bien la distribución del crecimiento, limitamos el mapa a un radio de 2,5 km desde el centro de la ciudad. Usamos la función geocode_OSM() del paquete tmaptools para obtener las coordenadas de Valencia en clase sf. Después proyectamos los puntos al sistema que usamos para el edificado (EPSG:25830). Como último paso creamos con la función st_buffer() un buffer con 2500 m y la intersección con nuestros datos de los edificios. También es posible crear un buffer en forma de un rectángulo indicando el tipo de estilo con el argumento endCapStyle = "SQUARE".

# obtenemos las coordinadas de Valencia
ciudad_point <- tmaptools::geocode_OSM("Valencia", 
                                      as.sf = TRUE)

# proyectamos los datos
ciudad_point <- st_transform(ciudad_point, 25830)

# creamos un buffer
point_bf <- st_buffer(ciudad_point, 2500)


# obtenemos la intersección entre el buffer y la edificación
buildings_val25 <- st_intersection(buildings_val, point_bf)

## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Preparar los datos para el mapas

Para poder visualizar bien las diferentes épocas de crecimiento, categorizamos el año en 15 grupos empleando cuartiles.

#encontrar 15 clases
br <- classIntervals(year(buildings_val25$beginning), 15, "quantile")

## Warning in classIntervals(year(buildings_val25$beginning), 15, "quantile"):
## var has missing values, omitted in finding classes

#crear etiquetas
lab <- names(print(br, under = "<", over = ">", cutlabels = FALSE))

## style: quantile
##      < 1890 1890 - 1912 1912 - 1925 1925 - 1930 1930 - 1940 1940 - 1950 
##         940        1369         971         596        1719        1080 
## 1950 - 1957 1957 - 1962 1962 - 1966 1966 - 1970 1970 - 1973 1973 - 1977 
##        1227        1266        1233        1165        1161         932 
## 1977 - 1987 1987 - 1999      > 1999 
##        1337        1197        1190

#categorizar el año
buildings_val25 <- mutate(buildings_val25, 
               yr_cl = cut(year(beginning), br$brks, labels = lab, include.lowest = TRUE))

Mapa de Valencia

El mapa creamos con el paquete tmap. Es una interesante alternativa a ggplot2. Se trata de un paquete de funciones especializadas en crear mapas temáticos. La filosofía del paquete sigue a la de ggplot2, creando multiples capas con diferentes funciones, que siempre empiezan con tm_* y se combinan con +. Cuando necesitamos más colores del máximo permitido por RColorBrewer podemos pasar los colores a la función colorRampPalette(). Esta función interpola para un mayor número más colores de la gama.

#colores
col_spec <- RColorBrewer::brewer.pal(11, "Spectral")

#función de una gama de colores
col_spec_fun <- colorRampPalette(col_spec)


#crear los mapas
tm_shape(buildings_val25) +
  tm_polygons("yr_cl", 
              border.col = "transparent",
              palette = col_spec_fun(15),
              textNA = "Sin dato",
              title = "") +
 tm_layout(bg.color = "black",
           outer.bg.color = "black",
           legend.outside = TRUE,
           legend.text.color = "white",
           legend.text.fontfamily = "Montserrat", 
            panel.label.fontfamily = "Montserrat",
            panel.label.color = "white",
            panel.label.bg.color = "black",
            panel.label.size = 5,
            panel.label.fontface = "bold")

Mapa dinámico en leaflet

Una ventaja muy interesante es la función tmap_leaflet() del paquete tmap para pasar de forma sencilla un mapa creado en el mismo marco a leaflet.

#mapa tmap de Santiago
m <-   tm_shape(buildings_val25) +
          tm_polygons("yr_cl", 
              border.col = "transparent",
              palette = col_spec_fun(15),
              textNA = "Sin dato",
              title = "")


#mapa dinámico
tmap_leaflet(m)